home *** CD-ROM | disk | FTP | other *** search
/ Mac-Source 1994 July / Mac-Source_July_1994.iso / C and C++ / Science⁄Math / VideoToolbox / VideoToolboxSources / ConvolveX.c < prev    next >
Text File  |  1993-07-09  |  11KB  |  285 lines
  1. /* 
  2. ConvolveX/Y.c
  3. © 1989-1993 Denis G. Pelli
  4. Simple one-dimensional convolution. Convolve image in srcBits by array f[] to
  5. produce image dstBits. src and dst may have any pixelSize, but the calculation
  6. itself is done with 8 bits per pixel. Assumes that everything has its origin at
  7. its center. Before calling ConvolveX/Y you should call SetEntries() to set up
  8. the color tables appropriately (unless source and destination pixmaps are both 8
  9. bits/pixel, in which case the color tables are ignored).
  10.  
  11. ConvolveX/Y fills each pixel in the destination with a weighted sum over the
  12. source, centered on the corresponding pixel in the source. To avoid edge
  13. effects, it only computes pixels in the destination which have full support in
  14. the source. If necessary, ConvolveX/Y will inset the source or destination rect,
  15. so that the destination width/height plus the array length equals the source
  16. width/height.
  17.  
  18. It would be nice to convert this to use Get and SetPixelsQuickly instead of
  19. directly accessing the pixels. This would introduce a small overhead, since the
  20. numbers would be copied into a long array and from thence into the pixels, instead
  21. of being stuffed directly, but would probably make this program quite a bit easier
  22. to read.
  23.  
  24. BUGS:
  25. Setting IMPROVED_ROUNDING to true improves results in some applications, but
  26. screws up others (e.g. the demo Filter.c). I don't have time to figure this out at 
  27. the moment, so I'm leaving a flag that can be turned on or off to suit your 
  28. application. I expect to fix this and eliminate the flag in the future.
  29.  
  30. HISTORY:
  31. 4/1/89 dgp    twice as fast, thanks to recoding to use long instead of float arithmetic.
  32. 4/2/89 dgp    made universal by using line buffers for both src and dst, and interfacing
  33.             to the src and dst PixMaps solely through CopyBits.
  34. 4/2/89 dgp    ConvolveX/Y.c, a pre-processor conditional was added to generate ConvolveX
  35.             and ConvolveY from the SAME source. The two files are IDENTICAL except
  36.             for the preprocessor definition of CONVOLVE_X as 1 to generate ConvolveX,
  37.             or 0 to generate ConvolveY.
  38. 10/9/89 dgp    Fixed overflow bug. Speeded up the convolution loop slightly.
  39.             Added conditional to use CopyBitsQuickly instead of CopyBits when source and
  40.             destination pixelSize are 8 bits. This has two virtues. It's faster, and
  41.             it bypasses the Color Manager's transformations through the color and inverse
  42.             color tables, which costs time and loses accuracy.
  43. 10/10/89 dgp Now use RectToAddress to determine the pixelSize of src and dst, to deal
  44.             correctly with window's pixmap.
  45. 7/25/91    dgp    Added explicit preprocessor symbol to turn diagnostic printout on and off.
  46. 8/24/91    dgp    Made compatible with THINK C 5.0.
  47. 8/27/92    dgp    Added Gestalt().
  48. 11/10/92 dgp Now initialize all fields of pixmap, especially the PMVersion.
  49.             Identify all error messages as originating from here: ConvolveX/Y.
  50.             Replaced compile-time conditional by a run-time conditional, so
  51.             that the ConvolveY.c file is no longer needed.
  52. 1/22/93    dgp    SwapMMUMode().
  53. 2/1/93    dgp Improved rounding. Removed calls to Gestalt, as unnecessary overhead.
  54. 2/7/93    dgp    Introduced switch IMPROVED_ROUNDING to disable rounding for debugging.
  55. 2/8/93    dgp Fixed overflow problem caused by the new rounding.
  56. 7/9/93    dgp check for 32-bit addressing capability.
  57. */
  58. #include "VideoToolbox.h"
  59. #include <math.h>
  60. void ConvolveXY(double f[],int dim,BitMap *srcBits
  61.     ,BitMap *dstBits,Rect *srcRectPtr,Rect *dstRectPtr,int convolveX);
  62. #define DIAGNOSTICS 0        // true or false
  63. #define IMPROVED_ROUNDING 1    // true or false
  64.  
  65. void ConvolveX(double f[],int dim,BitMap *srcBits
  66.     ,BitMap *dstBits,Rect *srcRectPtr,Rect *dstRectPtr)
  67. {
  68.     ConvolveXY(f,dim,srcBits,dstBits,srcRectPtr,dstRectPtr,1);
  69. }
  70.  
  71. void ConvolveY(double f[],int dim,BitMap *srcBits
  72.     ,BitMap *dstBits,Rect *srcRectPtr,Rect *dstRectPtr)
  73. {
  74.     ConvolveXY(f,dim,srcBits,dstBits,srcRectPtr,dstRectPtr,0);
  75. }
  76.  
  77. void ConvolveXY(double f[],int dim,BitMap *srcBits
  78.     ,BitMap *dstBits,Rect *srcRectPtr,Rect *dstRectPtr,int convolveX)
  79. {
  80.     Rect mysrcRect,mydstRect;
  81.     int dstLines,srcLines,dstLength,srcLength;
  82.     int lines;
  83.     unsigned char *dstPtr=NULL;
  84.     register unsigned char *srcPtr=NULL;
  85.     register int i;
  86.     Rect dstRect, srcRect;
  87.     register long t;
  88.     register long *ifun, *fPtr;
  89.     long ifunScale;
  90.     int x,y,srcx,imin,imax;
  91.     double tmax;
  92.     PixMap srcTmp,dstTmp;    /* line buffers */
  93.     #if DIAGNOSTICS
  94.         GDHandle device=NULL;
  95.         double a;
  96.     #endif
  97.     short srcPixelSize,dstPixelSize;
  98.     int error;
  99.     long value;
  100.     char mmuMode=true32b;
  101.     Boolean can32;
  102.     
  103.     value=0;
  104.     error=Gestalt(gestaltAddressingModeAttr,&value);
  105.     can32=!error && (value&(1<<gestalt32BitCapable));
  106.     RectToAddress((PixMap *)srcBits,srcRectPtr,NULL,&srcPixelSize,NULL);
  107.     RectToAddress((PixMap *)dstBits,dstRectPtr,NULL,&dstPixelSize,NULL);
  108.     ifun = (long *) malloc(dim*sizeof(*ifun));
  109.     if(ifun == NULL) {
  110.         PrintfExit("ConvolveX/Y: Sorry, couldn't allocate enough memory for working array\007\n");
  111.     }
  112.     /*
  113.     We will sum in a long register.
  114.     The sum over the point spread function f[] will yield at most tmax, defined below,
  115.     so we can safely multiply f[] by LONG_MAX/tmax, without danger of overflowing t.
  116.     */
  117.     for (tmax=0.0, i=0;i<dim;i++) tmax += fabs(f[i]);
  118.     tmax *= UCHAR_MAX+1;    /* maximum value for a pixel, 255, plus 1 for rounding */
  119.     if(tmax<1.0)tmax=1.0;
  120.     ifunScale = LONG_MAX/tmax;
  121.     for (i=0;i<dim;i++) ifun[i] = f[i]*ifunScale;
  122.     
  123.     #if DIAGNOSTICS
  124.         /* This is a diagnostic printout to track down overflow problems */
  125.         t = 0L;
  126.         a = 0.0;
  127.         for (i=0;i<dim;i++) {
  128.             t += labs(UCHAR_MAX*ifun[i]);
  129.             a += fabs(f[i]);
  130.         }
  131.         device=GetGDevice();
  132.         SetGDevice(GetMainDevice());
  133.         PmForeColor(UCHAR_MAX);
  134.         PmBackColor(0);
  135.         PrintfExit("a=%lf, t=%lx, ifunScale=%lx, t/ifunScale=%lf\n",a,t,ifunScale,t*1.0/ifunScale);
  136.         SetGDevice(device);
  137.     #endif
  138.     
  139.     /*
  140.     Unfortunately CopyBits will crash if the line buffer is too large, e.g.
  141.     640 pixels long. We can minimize this problem by reducing the size of the
  142.     srcRect to be no larger than necessary (i.e. dstRect plus dim).
  143.     */
  144.     mysrcRect = *srcRectPtr;
  145.     mydstRect = *dstRectPtr;
  146.     if(convolveX){
  147.         i=mysrcRect.right-mysrcRect.left - (dim+mydstRect.right-mydstRect.left);
  148.         if(i>0)InsetRect(&mysrcRect,i/2,0);
  149.         i=mysrcRect.bottom-mysrcRect.top - (mydstRect.bottom-mydstRect.top);
  150.         if(i>0)InsetRect(&mysrcRect,0,i/2);
  151.         if(i<0)InsetRect(&mydstRect,0,-i/2);
  152.     }else{
  153.         i=mysrcRect.right-mysrcRect.left - (mydstRect.right-mydstRect.left);
  154.         if(i>0)InsetRect(&mysrcRect,i/2,0);
  155.         if(i<0)InsetRect(&mydstRect,-i/2,0);
  156.         i=mysrcRect.bottom-mysrcRect.top - (dim+mydstRect.bottom-mydstRect.top);
  157.         if(i>0)InsetRect(&mysrcRect,0,i/2);
  158.     }
  159.     srcRectPtr = &mysrcRect;
  160.     dstRectPtr = &mydstRect;
  161.  
  162.     if(convolveX){
  163.         dstLines = dstRectPtr->bottom-dstRectPtr->top;
  164.         srcLines = srcRectPtr->bottom-srcRectPtr->top;
  165.         dstLength = dstRectPtr->right-dstRectPtr->left;
  166.         srcLength = srcRectPtr->right-srcRectPtr->left;
  167.     }else{
  168.         dstLength = dstRectPtr->bottom-dstRectPtr->top;
  169.         srcLength = srcRectPtr->bottom-srcRectPtr->top;
  170.         dstLines = dstRectPtr->right-dstRectPtr->left;
  171.         srcLines = srcRectPtr->right-srcRectPtr->left;
  172.     }
  173.     if(dstLines > srcLines) lines = srcLines; 
  174.     else lines = dstLines;
  175.     
  176.     /* Allocate our line buffers. */
  177.     srcTmp = **(*GetGDevice())->gdPMap;    /* Use current device to init PixMap fields */
  178.     srcTmp.pmVersion=srcTmp.packType=srcTmp.packSize=0;
  179.     srcTmp.planeBytes=srcTmp.pmReserved=0;
  180.     srcTmp.cmpCount=1;
  181.     srcTmp.pixelSize=8;                    /* The main loop assumes 8 bits per pixel */
  182.     if(convolveX){
  183.         SetRect(&srcTmp.bounds,0,0,srcLength,1);    /* a horizontal line buffer */
  184.     }else{
  185.         SetRect(&srcTmp.bounds,0,0,1,srcLength);    /* a vertical line buffer */
  186.     }
  187.     srcTmp.rowBytes = 2*(((srcTmp.bounds.right-srcTmp.bounds.left)*srcTmp.pixelSize+15)/16);
  188.     srcTmp.baseAddr = (Ptr) malloc(srcTmp.rowBytes*(srcTmp.bounds.bottom-srcTmp.bounds.top));
  189.     srcTmp.rowBytes |= 0x8000;    /* Mark it as a PixMap, not a BitMap */
  190.     dstTmp = srcTmp;
  191.     if(convolveX){
  192.         SetRect(&dstTmp.bounds,0,0,dstLength,1);    /* a horizontal line buffer */
  193.     }else{
  194.         SetRect(&dstTmp.bounds,0,0,1,dstLength);    /* a vertical line buffer */
  195.     }
  196.     dstTmp.rowBytes = 2*(((dstTmp.bounds.right-dstTmp.bounds.left)*dstTmp.pixelSize+15)/16);
  197.     dstTmp.baseAddr = (Ptr) malloc(dstTmp.rowBytes*(dstTmp.bounds.bottom-dstTmp.bounds.top));
  198.     dstTmp.rowBytes |= 0x8000;    /* Mark it as a PixMap, not a BitMap */
  199.     if(srcTmp.baseAddr == NULL || dstTmp.baseAddr == NULL)
  200.         PrintfExit("ConvolveX/Y: Sorry, couldn't allocate enough memory for working arrays\007\n");
  201.     
  202.     /* Define the rectangles that map src and dst to their respective line buffers */
  203.     srcRect = srcTmp.bounds;
  204.     OffsetRect(&srcRect,srcRectPtr->left,srcRectPtr->top);
  205.     dstRect = dstTmp.bounds;
  206.     OffsetRect(&dstRect,dstRectPtr->left,dstRectPtr->top);
  207.     if(convolveX){
  208.         OffsetRect(&srcRect,0,srcLines/2-lines/2);
  209.         OffsetRect(&dstRect,0,dstLines/2-lines/2);
  210.     }else{
  211.         OffsetRect(&srcRect,srcLines/2-lines/2,0);
  212.         OffsetRect(&dstRect,dstLines/2-lines/2,0);
  213.     }
  214.     
  215.     /* Convolve one line per iteration */
  216.     /* y indicates which line, x indicates which pixel along that line */
  217.     for (y= -lines/2; y<(lines+1)/2; y++) {
  218.         /* first copy a line from srcBits into line-buffer srcTmp */
  219.         if(srcPixelSize == srcTmp.pixelSize)
  220.             CopyBitsQuickly(srcBits,(BitMap *) &srcTmp,&srcRect,&srcTmp.bounds,srcCopy,NULL);
  221.         else
  222.             CopyBits(srcBits,(BitMap *) &srcTmp,&srcRect,&srcTmp.bounds,srcCopy,NULL);
  223.         if(can32)SwapMMUMode(&mmuMode);
  224.         if(!convolveX){
  225.             /* Unfortunately the minimum PixMap rowBytes is 2 bytes, so we
  226.             squeeze out the extra space before processing
  227.             */
  228.             srcPtr= (unsigned char *) srcTmp.baseAddr;
  229.             dstPtr= (unsigned char *) srcTmp.baseAddr;
  230.             for (x=0; x<srcLength; x++) {
  231.                 *dstPtr++ = *srcPtr;
  232.                 srcPtr += 2;
  233.             }
  234.         }
  235.         dstPtr= (unsigned char *) dstTmp.baseAddr;
  236.         srcx = srcLength/2 - dim/2 - dstLength/2;
  237.         for (x=0; x<dstLength; x++) {
  238.             imax=dim;
  239.             if(srcx+dim>srcLength) imax=srcLength-srcx;
  240.             if(srcx<0){
  241.                 imin= -srcx;
  242.                 srcPtr=(unsigned char *) srcTmp.baseAddr;
  243.             }
  244.             else {
  245.                 imin=0;
  246.                 srcPtr=(unsigned char *) srcTmp.baseAddr+srcx;
  247.             }
  248.             fPtr=&ifun[imin];
  249.             for (t=0L, i=imin;i<imax;i++) t += *fPtr++ * (long)(unsigned long) *srcPtr++;
  250. #if IMPROVED_ROUNDING
  251.             if(t>0)t+=ifunScale>>1;        // round t/ifunScale to nearest integer
  252.             else t-=ifunScale>>1;
  253. #endif
  254.             t/=ifunScale;
  255.             if(convolveX){
  256.                 *dstPtr++ = (unsigned char) t;
  257.             }else{
  258.                 /* Unfortunately the minimum PixMap rowBytes is 2 bytes, so we
  259.                 insert a space between pixels
  260.                 */
  261.                 *dstPtr = (unsigned char) t;
  262.                 dstPtr += 2;
  263.             }
  264.             srcx++;
  265.         }
  266.         if(can32)SwapMMUMode(&mmuMode);
  267.         /* finally, copy the line from buffer dstTmp to dstBits */
  268.         if(dstPixelSize == dstTmp.pixelSize)
  269.             CopyBitsQuickly((BitMap *) &dstTmp,dstBits,&dstTmp.bounds,&dstRect,srcCopy,NULL);
  270.         else
  271.             CopyBits((BitMap *) &dstTmp,dstBits,&dstTmp.bounds,&dstRect,srcCopy,NULL);
  272.         if(convolveX){
  273.             OffsetRect(&srcRect,0,1);    /* adjust Rects to point to next line */
  274.             OffsetRect(&dstRect,0,1);
  275.         }else{
  276.             OffsetRect(&srcRect,1,0);    /* adjust Rects to point to next line */
  277.             OffsetRect(&dstRect,1,0);
  278.         }
  279.     }
  280.     free((void *) ifun);
  281.     free((void *) srcTmp.baseAddr);
  282.     free((void *) dstTmp.baseAddr);
  283. }
  284.  
  285.